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Abstract 

The discovery of the chaotic behavior of the planetary 
orbits in the Solar System (Laskar, 1989, 1990) was ob- 
tained using numerical integration of averaged equations. 
In (Laskar, 1994), these same equations arc integrated 
over several Gyr and show the evidence of very large possi- 
ble increase of the eccentricity of Mercury through chaotic 
diffusion. On the other hand, in the direct numerical in- 
tegration of (Ito and Tanikawa, 2002) performed without 
general relativity over ±4 Gyr, the eccentricity of Mer- 
cury presented some chaotic diffusion, but with a max- 
imal excursion smaller than about 0.35. In the present 
work, a statistical analysis is performed over more than 
1001 different integrations of the secular equations over 5 
Gyr. This allows to obtain for each planet, the probabil- 
ity for the eccentricity to reach large values. In particular, 
we obtain that the probability of the eccentricity of Mer- 
cury to increase beyond 0.6 in 5 Gyr is about 1 to 2 %, 
which is relatively large. In order to compare with (Ito 
and Tanikawa, 2002), we have performed the same anal- 
ysis without general relativity, and obtained even more 
orbits of large eccentricity for Mercury. In order to clar- 
ify these differences, we have performed as well a direct 
integration of the planetary orbits, without averaging, for 
a dynamical model that do not include the Moon or gen- 
eral relativity with 10 very close initial conditions over 
3 Gyr. The statistics obtained with this reduced set are 
comparable to the statistics of the secular equations, and 
in particular we obtain two trajectories for which the ec- 
centricity of Mercury increases beyond 0.8 in less than 1.3 
Gyr and 2.8 Gyr respectively. These strong instabilities 
in the orbital motion of Mecury results from secular reso- 
nance beween the perihelion of Jupiter and Mercury that 
are facilitated by the absence of general relativity. The 
statistical analysis of the 1001 orbits of the secular equa- 
tions also provides probability density functions (PDF) 
for the eccentricity and inclination of the terrestrial plan- 
ets (Mercury, Venus, the Earth and Mars) that are very 



well approximated by Rice PDF. This provides a very sim- 
ple representation of the planetary PDF over 5 Gyr. On 
this time-scale the evolution of the PDF of the terrestrial 
planets is found to be similar to the one of a diffusive 
process. As shown in (Laskar, 1994), the outer planets 
orbital elements do not present significant diffusion, and 
the PDFs of their eccentricities and inclinations are well 
represented by the PDF of quasiperiodic motions with a 
few periodic terms. 
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1 Introduction 

The discovery of the chaotic behavior of the planetary or- 
bits in the Solar System (Laskar, 1989, 1990) was obtained 
using numerical integration of averaged equations. Since 
then, this chaotic behavior has been confirmed through 
direct numerical simulation, without averaging (Quinn et 
al, 1991, Sussman and Wisdom, 1992). More recently, 
integrations of more accurate planetary equations have 
been performed over 100 to 250 Myr (Varadi et ai, 2003, 
Laskar et ai, 2004a, b). Due to the chaotic evolution 
of the system, and to the uncertainty on the model and 
initial conditions, the interval of precise validity of these 
solutions is limited to about 40 Myr (Laskar et al, 2004a), 
even if some components of the solutions can be used over 
longer time for paleoclimate studies, as the 405 kyr oscil- 
lation of the Earth eccentricity (Laskar et al, 2004a). 

Over periods of time longer than 100 Myr, it becomes 
hopeless to search for a precise solution for the orbital pa- 
rameters of the inner planets (Mercury, Venus, Earth and 
Mars). On the other hand, it is important to understand 
the possible behavior of these solutions, and in particu- 
lar of the possible variations of the action variables of the 
orbits (semi major axis, eccentricity and inclination). 

1 E-mail address: laskar@imcce.fr 



A first study of the chaotic diffusion of the planetary 
orbits was made in (Laskar, 1994) is order to search for 
the maximum possible variations of the eccentricities and 
inclinations of the planets, using the secular equations. 

In the present work, a statistical view of the chaotic 
evolution of the planetary orbits is described for all plan- 
ets of the Solar System (Pluto is excluded) . This study is 
based on the numerical integrations of 1001 different solu- 
tions of the averaged equations of the Solar System using 
very close initial conditions, compatible with our present 
knowledge. Some results concerning the orbital evolution 
of Mars were already presented in (Laskar et ai, 2004b) 
and we will denote this paper PI in the following, since in 
many cases, we will refer to this previous paper to avoid 
duplication. 

2 The secular equations 

In order to investigate the diffusion of the orbits over 5 
Gyr, we will use the secular equations of Laskar (1990), 
with some small modifications. The secular equations are 
obtained by averaging the equations of motion over the 
fast angles that are the mean longitudes of the planets. 
The averaging of the equation of motion is obtained by 
expanding the perturbations of the Keplerian orbits in 
Fourier series of the angles, where the coefficients them- 
selves are expanded in series of the eccentricities and in- 
clinations. This averaging process was conducted in a 
very extensive way, up to second order with respect to 
the masses, and through degree 5 in eccentricity and in- 
clination, leading to truncated secular equations of the 
Solar System of the form 

= V^l{Tw + <I> 3 (w,w) + <f> 5 (w,w)} (1) 

where w = (z 1 , . . . ,z 8 ,Ci, ■ ■ • ,Cs), with z k = e k exp(tu fe ), 
Qu = sin(«fc/2) exp(fifc) (vjk is the longitude of the perihe- 
lion). The 16 x 16 matrix T is the linear Lagrange-Laplace 
system, while ^(iu, w) and <&5(w,w) gather the terms of 
degree 3 and 5. 

The system of equations thus obtained contains some 
150000 terms, but can be considered as a simplified sys- 
tem, as its main frequencies are now the precession fre- 
quencies of the orbits of the planets, and no longer com- 
prise their orbital periods. The full system can thus be 
numerically integrated with a very large step-size of 200 
to 500 years. Contributions due to the secular perturba- 
tion of the Moon and general relativity are also included 
(see Laskar 1990, 1996, and laskar et ai, 2004b (PI) for 
more details and references). 

This secular system is then simplified and reduced to 
about 50000 terms, after neglecting terms of very small 



value (Laskar 1994). Finally, a small correction of the 
terms of the matrix T of ([!]), after diagonalization, is per- 
formed in order to adjust the linear frequencies, in a sim- 
ilar way as it was done in (Laskar 1990). Indeed, in the 
outer planetary system, terms of higher order are of some 
importance, but their main effect will be to slightly mod- 
ify the values of the main frequencies of the system. The 
correction that is done here is a simple way to correct 
for this effect. With the present small adjustment, the 
secular solutions are very close to the direct numerical in- 
tegration La2004 (Laskar et ai, 2004a, b) over about 35 
Myr (PI, Figs. 15 and 16). As noted in PI, this time is 
about the time over which the direct numerical solution 
itself is valid (Laskar et ai, 2004a, Figs. 20, 21), be- 
cause of the imperfections of the model. Moreover, as the 
step-size used in the secular equations is 200 years instead 
of 1.82625 days for the direct integration, over very long 
times the numerical noise will be smaller. It is thus le- 
gitimate to investigate the diffusion of the orbital motion 
over long times using the secular equations. The major 
advantage, besides reducing the roundoff errors, resides 
in the computation speed : the integration of the sec- 
ular equations is 2000 times faster than the integration 
of the non-averaged equations, and we can compute a 5 
Gyr solution for the Solar System in 12 hours on a Com- 
paq alpha workstation (833 Mhz). We are thus able to 
make statistics over many solutions with close initial con- 
ditions. In these computations, our main limitation will 
be the huge amount of data generated by these numerical 
integrations. 
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Table 1: Offsets of the initial conditions for the 1001 integra- 
tions of the secular system for the variables k = ecosro and 
h = esinro. The different integrations correspond to offsets of Ne 
for N = —100, . . . ,+100 and e = 10~ 10 , in a single variable, for 
a single planet, while the other variables are kept to their nominal 
values. 



3 Maximum excursion 

For the present analysis, we have integrated 1001 orbital 
trajectories of the Solar System over 5 Gyr in negative 
time with very small variations of the initial conditions 
with respect to the nominal solution. The initial condi- 
tions of the secular system nominal solution are the same 
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Figure 1: Probability (in percents) to reach a given value of the eccentricity (left) or inclination (right) over a given time interval. The 
statistic is established with 1001 solution with very close initial conditions. The maximum values are computed over 50, 100, 250, 500, 
1000, 2000, 3000, 4000, and 5000 Myr for the inner planets. For the outer planets, the diffusion is so small that these computations are 
easily summarized in a Table. Inclinations are computed with respect to the ecliptic J2000. 



as Table 1 and 2 from (Laskar, 1986), and are derived from 
the initial conditions of the VSOP82 solution of (Bre- 
tagnon, 1982). The phase space of the secular system 
is of real dimension 32. The various initial conditions 
for the 1001 cases are obtained with a small variation of 
the initial value of a single secular variable k = ecos(ro) 
or h — esin(n7), for a single planet, according to Table [I] 
leaving the other 31 initial variables equal to their nomi- 
nal values. 

As we have performed 1001 numerical simulations of 
the whole Solar System over 5 Gyr, it is impossible to 
display the detailed results of these integrations, and we 
have chosen here to describe the most significant features 
of these integrations. One important point, that was ad- 



dressed in (Laskar, 1994) is the maximum value reached 
by the eccentricity or inclination, as a result of the chaotic 
diffusion of the trajectories. In (Laskar, 1994), only 5 so- 
lutions were followed, after some small change of initial 
conditions in every interval of 500 Myr. This was an eco- 
nomical way to obtain the maximal possible value for each 
parameter, but without any estimate of the probability 
to reach these values. In particular, I could demonstrate 
that it was possible for Mercury to reach very high values 
for its eccentricity, allowing eventually a close encounter 
with Venus, while similar crossings were not reached for 
the other planets. 

Here we perform more conventional statistics, and all 
1001 integrations have been followed over the whole 5 
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Table 2: In negative time, number of integrated solutions (out of a 
total of 1001 cases, rescaled to 1000) for which the maximum value 
reached by the eccentricity of Mercury (e max ) over a given time 
(500Myr, 1000 Myr, 1500 Myr, 2000 Myr, 3000 Myr, 4000 Myr, or 
5000 Myr) is above a specified value (e m o). 
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Table 3: For positive time , number of integrated solutions (out of 
a total of 478 cases, rescaled to 1000) for which the maximum value 
reached by the eccentricity of Mercury (e max ) over a given time 
(500Myr, 1000 Myr, 1500 Myr, 2000 Myr, 3000 Myr, 4000 Myr, or 
5000 Myr) is above a specified value (e m n)- 

Gyr interval. As it was mentioned before, (Laskar, 1990, 
1994), there is practically no diffusion for the outer planet 
system that behaves nearly as a quasiperiodic and regular 
system. On the opposite, there is a significant diffusion 
of the eccentricities and inclinations of the inner planets. 
The statistics on the maximum values reached by the ec- 
centricity and inclination over different time intervals of 
50, 100, 250, 500, 1000, 2000, 3000, 4000, and 5000 Myr 
are displayed in figure [l] for Mercury, Venus, The Earth 
and Mars. 

Over 50 Myr, the effect of the diffusion is not yet no- 
ticeable, and all solutions reach the same maximum value. 
This is reflected by a vertical curve in figure [T] But be- 
yond this time interval, significant differences appear. In 
order to make these statistics more readable, the lower 
part of these maximum graphs, corresponding to the 5 
percents of the solutions with largest variations of the ec- 
centricities and inclination has been enlarged 

3.1 Discussion 

An important aspect of this study is the estimate of the 
probability for the eccentricity of Mercury to reach high 
values, allowing a possible close encounter with Venus. 



Over 1 Gyr, all solutions remained with e max < 0.5, and 
only 0.1% of the solutions went beyond 0.5 over 2 Gyr 
(Table |2|. On the other hand, over longer time inter- 
val, the chaotic diffusion allowed Mercury's eccentricity 
to reach very high values. In our simulations, 0.9% of the 
solutions reached an eccentricity of 0.9 within 5 Gyr, and 
0.6% within 4 Gyr. At this point, one should remind that 
a close encounter of Mercury and Venus is only possible 
if the eccentricity of Mercury reaches values of about 0.75 
(assuming an eccentricity of Venus of 0.06). We still have 
in our simulations about 1 % of the solutions that would 
allow for a close encounter with Venus. 

In order to test the stability of the results displayed in 
Table [2j I have also plotted in Table [3] the results obtained 
with a similar experiment, but performed in positive time. 
The equations are the same, but the integration step is 
now positive. The simulation is made over 478 different 
orbit^] but the results are scaled to 1000 for easier com- 
parison with Table [2] In a same way as for a diffusive 
process, the results in positive time are very similar to 
the results obtained in negative time, which is what was 
expected. 

In doing these estimates, we need to keep in mind 
that the equations that are integrated here are the av- 
eraged equations of motions, where the disturbing func- 
tion of the mutual interactions of the planets is expanded 
in series of eccentricities and inclinations. This expan- 
sion is divergent for high eccentricities. Indeed, already 
in the two-body problem, the expansion of the eccentric 
anomaly in powers of excentricity has a radius of conver- 
gence e m « 0.6627 (see Wintner, 1947). Additionnally, 
the inverse of the distance becomes singular at collision, 
that is for an eccentricity of Mercury of about 0.75, de- 
pending on the eccentricity of Venus. Moreover, in the 
vicinity of collisions, numerous mean motion resonances 
overlap, giving rise to mean motion chaotic behavior, with 
possible changes in the semi-major axis of the planets, 
while these are constant in the secular system (see Fig. 

Nevertheless, although the truncations of series expan- 
sions involved in the construction of the secular system 
are made without estimates of the remainders, as it is 
usually the case in astronomy, I conjecture here that up 
to an eccentricity of Mercury of about 0.6, the dynamics 
of the full system is well represented by the dynamics of 
the secular system. Indeed, it has been observed that in 
many cases, the range of validity of secular systems ex- 
tends much beyond theoretical estimates (see for example 

x This odd number of cases is purely accidental. Our dedicated 
parallel computer broke down, and as this positive time integrations 
are used only as a check, we considered that 478 cases were sufficient 
for this purpose. 
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Libert and Henrard, 2006). 

Moreover, my assumption is that in general, the full non 
averaged system is less regular than the secular system. I 
can thus assume that for eccentricities of Mercury below 
0.6 the data displayed in Table |2]provide good estimates of 
the actual probabilities to reach these eccentricity values, 
while for eccentricities above 0.6, the data given in Table 
[2] are only lower bounds of these probabilities. 

It is certainly desirable that the same statistical studies 
should be conducted with the full equations of motions, 
without averaging, although this will require a consider- 
able amount of CPU time that is still difficult to obtain. 

On the other hand, the data displayed in Table [2] 
provide substantially more information than in (Laskar, 
1994) where only the possibility of reaching high values of 
Mercury's eccentricity allowing for a collision with Venus 
was demonstrated. Here we show that the probability to 
reach these high values is relatively large (about 1 to 2%). 
It thus becomes conceivable to prepare in the near future 
a full scale numerical simulation of the same problem with 
the non averaged equations. 

3.2 Comparison with direct integration 
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Figure 2: Maximum value reached by the relative variation of the 
total angular momentum (|(c(i) — c(0))/c(0)|) of the Solar System 
for the 10 solutions S\, .. -Siq. The maximum values (dcmax) are 
computed over 10 Myr intervals. 

It would be interesting to compare with an integration 
of the full equations of motion over Gyr time scale, but 
the author is not aware of a single numerical integration of 
the kind with a comparable dynamical model, including 
the contribution of the Moon and of general relativity. 
The long term integrations of (Varadi et al. 2003, Laskar 
et ai, 2004a,b) use precise models but span only a few 
100 Myr. On the other hand, the long term integration 
of Ito and Tanikawa (2002) is performed over a few Gyr 
but do not include the relativistic contribution. 
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Figure 3: Maximum value reached by the eccentricity of Mercury 
over 3Gyr. The maximum values (emax) are computed over 500 
Myr intervals for the 10 solutions Si,.. . Sio- 
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Figure 4: Maximum value reached by the relative variation of the 
total energy (\{h(t) - h(0))/h(Q)\) of the Solar System for the 10 
solutions Si, . . . iSio- The maximum values (dhmax) are computed 
over 10 Myr intervals. 



3.2.1 The integration of Ito and Tanikawa 

The numerical integration of Ito and Tanikawa (2002) 
is an integration of the Newtonian planetary equations, 
without relativistic contributions. It does not comprise ei- 
ther the effect of the Moon as a separate body. In fact, the 
mass of the Moon has been added to the mass of the Earth 
for five integrations (JV+i, 7V+2, N-i,N-2) that span 
from 3.9 to 5 Gyr in the past (iV_i, N-2) or in the future 
(iV-|_i, iV+2, iV+3). The initial conditions are taken from 
the JPL ephemeris DE245 (Standish, 1990). The integra- 
tor is a second order symplectic integrator (Wisdom and 
Holman 1991) with a step size of 8 days. From the figure 
1 of (Ito and Tanikawa 2002) we can estimate that the 
maximum relative error in angular momentum is about 
5 x 10~ n for all solutions except for one for which it is 
about 20 x 10 -11 , the difference being probably due to dif- 
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ferent hardware. The maximum relative error in energy 
is about 6 x 10~ 9 . 

The authors do not provide very precise details on the 
maximum values reached by the eccentricities, but men- 
tion that the general behavior of the eccentricity of Mer- 
cury is similar with the results of (Laskar, 1994, 1996), 
although they obtained a maximum eccentricity for Mer- 
cury of about 0.35 over ±4 Gyr. 

If we consider the numerical experiment that we con- 
ducted here with the secular system in negative time (Ta- 
ble [2| over 1001 solutions, we have only 32.8 % of the 
solutions that lead to an eccentricity of Mercury larger 
than 0.35, while 71.4 % had a maximal eccentricity larger 
than 0.32, and only 26.4 % an eccentricity larger than 
0.36. There is indeed a very rapid decrease of the proba- 
bility to reach a given eccentricity in the vicinity of 0.35. 

Due to the lack of precision on the maximal value 
reached by the eccentricity of Mercury in (Ito and 
Tanikawa 2002), we can consider that the event (e max <~ 
0.35) reached in their simulations has, according to our 
experiment on the secular system, a probability of about 
75% to occur. As they made 5 simulations, the resulting 
probability would be (0.75) 5 , that is about 24%. This is 
not very large, but not unrealistic. 

Nevertheless, there is a difference in the two experi- 
ments, as our secular equations comprise the relativistic 
contributions. One should thus also test the behavior of 
the secular system in absence of the relativistic contribu- 
tion. 

3.2.2 The secular system without relativity 

I have thus repeated the previous simulations in absence 
of relativity, expecting to find a more stable system. But 
the result was the opposite, and most of the solutions 
lead to high values of the eccentricity in 5 Gyr (Table 
[4|. One can see that in this case, the probability of the 
event (e max <~ 0.35) becomes less than 25 %, and the 
probability of having 5 solutions with this behavior (if 
they are not related) becomes totally unrealistic with a 
value smaller than (0.25) 5 , that is about 0.1%. 

3.3 A new direct integration without rel- 
ativity 

In order to clarify the situation, it was thus necessary 
to make a new direct numerical integration. Indeed, the 
probability of reaching high values of the eccentricity of 
Table [4] are so high that one should be able to obtain 
solutions with high eccentricity with a moderate number 
or trials. I thus decided to integrate 10 orbits, Si,..., Sio, 
with very close initial conditions over 3 Gyr. 
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Table 4: When General Relativity is not taken into account, num- 
ber of integrated solutions in negative time (out of a total of 988 
cases, rescaled to 1000) for which the maximum value reached by 
the eccentricity of Mercury (e m ax) over a given time (500Myr, 1000 
Myr, 1500 Myr, 2000 Myr, 3000 Myr, 4000 Myr, or 5000 Myr) is 
above a specified value (e m o). 
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Table 5: For the direct integration of Newtonian equations, . number 
of integrated solutions in negative time, (out of a total of 10 cases) 
for which the maximum value reached by the eccentricity of Mercury 
(e ma x) over a given time (500Myr, 1000 Myr, 1500 Myr, 2000 Myr, 
3000 Myr) is above a specified value (e m o)- 

The model comprises the 8 major planets and the dwarf 
planet Pluto with Newtonian interactions. The Earth- 
Moon barycenter with the sum of the masses of the Earth 
and Moon is used instead of the Earth. The integrator is 
the SABA4 symplectic integrator of (Laskar and Robutel, 
2001) that was already used in (Laskar et al., 2004a, b) 
with a step size of 2.5 x 10~ 2 years. For a perturbed 
Hamiltonian of the form H = A+eB, using this integrator 
with a step r is equivalent to integrate exactly a close by 
Hamiltonian H where the error of method H — H is of the 
order 0(r 8 e) + 0(r 2 e 2 ). The integration is conducted in 
extended precision on Itanium II processors with 80 bit 
arithmetics. 

All parameters and initial conditions of the nominal so- 
lution are the same as the ones used in the new high pre- 
cision planetary ephemeris INPOP06 that has been de- 
velopped in our group. The reader should refer to the 
associated publication (Fienga et al, 2008) and web site 
(www.imcce.fr/inpop) for more details. The 10 different 
solutions Sk (fc = 1, • • ■ , 10) have the same initial condi- 
tions as the nominal solution, except for a shift of k x 10 
cm in the initial position of the Earth. 

The maximum values reached by the eccentricity of 
Mercury in these 10 solutions have been plotted in figure 
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Figure 5: Example of very high value reached by the exentricity of 
Mercury in the S$ solution. The eccentricity is given versus time 
in Myrs (top). In the bottom frame, the relative variation of total 
energy (dh = \ (h(t) — h(0))/h(0)\) of the system is given. The step 
size of the integration is 1 X 10 — 2 yr on the interval [—1.284, —1.215] 
Gyr and 1 X 10 -3 yr over the interval [-1.2868, -1.284] Gyr. 



[3] One can see that the present results differ substan- 
tially from the results of Ito and Tanikawa (2002), as the 
eccentricity can reach much higher values than in their 
simulations. In particular, for 5 solutions, the eccentric- 
ity went beyond 0.4, and for two of them, 52 and S§, the 
eccentricity increased to very high values, beyond 0.8, in 
respectively 2776 and 1286.8 Myr (Fig. |J. 

In the symplectic integration, the total angular momen- 
tum is conserved. The relative variation of the total an- 
gular momentum is thus an estimate of the roundoff error 
in the integration. For all solutions, this error remains be- 
low 10 -13 over the total length of the integration (Fig{2j). 
The relative error in energy is in general below 10 , 
except when the eccentricity of Mercury reaches high val- 
ues (FigQ. For the solution S 5 that reaches very high 
values of the eccentricity, the step size was changed to 
1 x 10~ 2 yr on the interval [-1.284, -1.215] Gyr, and then 
to 1 x 10~ 3 yr over the interval [-1.2868, -1.284] Gyr in 
order to achieve a good conservation of the total energy. 
With these settings, even for an eccentricity of 0.8, the 
relative error of the energy remains below 3 x 10 -10 over 



the whole interval of integration (Fig. [5j. 

When we gather the results of these 10 numerical simu- 
lations in a table (Tab. [5]) similar to the tables [2j [3] [4] one 
can see that in our numerical results, the number of high 
eccentricities are somewhat lower, but quite comparable 
to the results of our table [4] especially if we consider the 
small number of events (10) on which the present statis- 
tics are made. 

The present results, that are conducted with slightly 
more accurate integrations than the ones of Ito and 
Tanikawa (2002), demonstrate clearly that the large 
excursions of the eccentricity of Mercury described in 
(Laskar, 1994) are actually present in the full integra- 
tion of the Newtonian equations. Moreover, the statistics 
obtained with the secular equations in absence of general 
relativity (Tab. [4} are very close to the same statistics 
obtained with the direct integration (Tab. |5j. 

We can thus assume that the statistics obtained over 
1001 integrations of the secular equations (Tab. [2| are 
representative of the full system. Actually, as was al- 
ready stated in (Laskar, 1994), we expect even that the 
full equations of motion will be slightly more unstable 
than the secular equations, especially for high eccentric- 
ities, when the overlap of mean motion resonances will 
increase the chaoticity of the orbits. 

3.4 The effect of relativity 

The difference of behavior of the secular system with (Ta- 
ble^ and without gener al relativity (GR) (Tabled is im- 
pressive. The results of the direct integrations performed 
without GR are also in good agreement with the results of 
the secular equation in absence of GR. The contribution 
of GR is thus essential in order to ensure the relative sta- 
bility of Mercury. It is therefore important to understand 
the contribution of GR in this problem. 

The most obvious effect of GR is to increase the perihe- 
lion frequencies of the planets and especially of Mercury. 
Indeed, the value of the secular frequency g\ related to 
Mercury is 5.59"/yr in the vicinity of the origin with a 
contribution of 0.43"/yr from general relativity. The con- 
tribution of GR to the perihelion velocity decreases to 
0.086"/yr for Venus, 0.038"/yr for the Earth, 0.013" /yr 
for Mars, and only 0.0006"/yr for Jupiter (e. g. Laskar, 
1999, Table 4). The main frequency of the longitude of 
perihelion of Jupiter, g-, ~ 4.257"/yr (Laskar et al, 2004a, 
Table 3) is thus not changed much by GR and the main 
effect of GR is to increase the difference g\ — g$ from 
0.90"/yr in absence of GR to 1.33"/yr in presence of GR. 

In order to analyse the effect of this change in the dy- 
namics of Mercury, I have performed several integrations 
of the secular system, with and without GR, for various 
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Figure 6: Secular resonance g± = 35. For initial eccentricity of 
Mercury (eO) from to 0.95 with step size 0.001, a numerical inte- 
gration of the secular system is performed over 40 Myr. The value 
of gi is obtained by frequency analysis over 20 Myr intervals and 
is plotted versus initial eccentricity eO in presence (al) or absence 
(a2) of GR. The horizontal line corresponds to the secular frequency 
35 ~ 4.257"/yr. In (bl) and (b2) are reported the maximum values 
of the eccentricity of Mercury reached over the 40 Myr interval, with 
(bl) and without (b2) relativity. With GR, gi is far from the secu- 
lar resonance gi = #5 and the variation of eccentricity is moderate. 
On the opposite, when GR is not considered, the secular frequency 
gi is smaller, and as the eccentricity of Mercury increases, the orbit 
gets trapped into the secular resonance g± = g§ that may drive the 
eccentricity to very high values (b2). 



values of the initial eccentricity of Mercury. The initial 
eccentricity of Mercury in the nominal solution is 0.2056, 
while in the present experiment, this value is varied with 
a step of 0.001 from to about 0.95 until the integra- 
tion crashes rapidly. The integrations are computed over 
40 Myr. A frequency analysis is performed in order to 
compute the secular frequencies gk , Sfc of the system. For 
each orbit, This frequency analysis (Laskar, 1990, 1993) 
is made over slidings 20 Myr intervals with a 1 Myr offset. 
The values obtained for g\ are reported in figure [6] with 
(al) and without (a2) GR. 

The differences between the two plots are striking. 
With GR, when the eccentricity increases, g\ is modi- 
fied, but mostly increases and although gi presents large 
variations resulting from the chaotic behavior of the sys- 
tem, the system never gets trapped into the g\ = g$ res- 
onance. On the opposite, in absence of GR, the g\ val- 
ues are smaller, and as the eccentricity increases beyond 
eO w 0.6, there exists a large chaotic zone related to the 
secular resonance g\ — g§, the g§ secular frequency be- 
ing plotted as an horizontal line in Figs [6^il,a2. When 
the system is in this resonance, the eccentricity of Mer- 
cury is driven to very large excentricities, beyond 1 (the 
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Figure 7: Evolution of the difference of longitude of perihelion (roi — 
■cos) °f Mercury and Jupiter for the 1S5 solution. The angle is plotted 
modulo 67r in order to better visualize the trapping into resonance 
that occurs from —1284 Myr to —1286 Myr, and corresponds to the 
strong increase of the eccentricity of Mercury from 0.4 to 0.8 (Fig. 



eccentricity is not bounded by 1 in the secular system). 
In figure [6)32, one trajectory starting at about 0.5 also in- 
creases to high values. Indeed, most probably, the chaotic 
region extends in the 0.4 — 0.5 region of eccenricity, but it 
will take some additional time to reach the region of very 
strong chaos beyond e = 0.6. 

This study of the secular system is made here to explore 
rapidly the phase space of the system, and to guide our 
intuition. Once we see the importance of the g\ — g§ res- 
onance in the secular system, we can test its contribution 
on the full equations with minimal computations. 

We can now verify that the large increase of eccentricity 
of the solution <S 5 plotted in figure [5] is actually triggered 
by the secular resonance g\ = 55. In figure [7J is plot- 
ted the difference W\ — of the longitudes of perihelion 
of Mercury {wij and Jupiter (n7 5 ) from —1260 Myr to 
— 1286.8 Myr. It appears that this angle is circulating un- 
til the two perihelions get locked from about —1284 Myr 
to —1286 Myr, that is over the time interval that corre- 
sponds to a steady increase of Mercury's eccentricity. In 
a similar way, the solution 5 2 presents as well a large in- 
crease of Mercury's eccentricity from 0.5 to about 0.9 in 
only 2.5 Myr that corresponds to a locking of the peri- 
helions of Mercury and Jupiter from —2773.5 to —2776 
Myr. 

It should be noted that for Mercury and Jupiter, and 
contrarily to the Earth, the relation with the longitudes 
of perihelion vj\ , 1375 and secular frequencies <?i , 55 , respec- 
tively, is straightforward. Indeed, in both cases, the secu- 
lar frequency gk is clearly the leading periodic term of the 
quasiperiodic expansion of Zk = ek exp(itUk) (see Table II 
from Laskar, 1990). 

In order to see the real effectiveness of this gi = g& 
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Figure 8: Evolution of the eccentricity (top) and semi-major axis 
(bottom) of Mercury in a numerical integration of the full Solar Sys- 
tem with the non averaged equations with general relativity where 
the PPN parameters 27-/3 has been artificially set to —2 instead 
of 1 in the standard setting. With this setting, and the present ini- 
tial conditions and parameters from INPOP06 (Fienga et al., 2008), 
the system is in the gi = 35 secular resonance. The eccentricity of 
Mercury then increases to very high values in only 4 Myr, where it 
reaches the region of mean motion resonances overlap where chaotic 
variations of the semi-major axis occur. 



resonance, I have performed an additional single numeri- 
cal experiment, with the full equations of motion, includ- 
ing GR, and starting with the nominal intitial conditions, 
that is the initial conditions of the planetary cphcmcris 
INPOP06 (Fienga et al, 2008), without any change of 
the initial eccentricity of Mercury. But in order to set 
the system into the g\ — g$ resonance from the starting 
time of integration, I have changed the value of the post 
Newtonian parameter f3 to 3 instead of (3 — 1 in standard 
GR, while keeping 7 = 1 (Will, 2006). The contribution 
factor of perihelion shift 27 — (3 is then —2 instead of +1. 
The effect of this modified relativity is thus to artificially 
decrease the secular frequency g\ by 0.86"/yr, instead of 
adding 0.43"/yr as in standard GR. The system is thus 
from the begining in the g\ = g 5 resonance, and the effect 
is immediate as Mercury's eccentricity increases beyond 
0.7 in less than 4 Myr (Fig. [8]) . Once in this region of high 
eccentricity, strong chaotic behavior due to short period 
resonances induces significant changes in the semi-major 
axis of Mercury (Fig. [8]). 



4 Density functions 

The maximal possible value of the planet eccentricity is 
important for the analysis of the system stability, but it 
concerns the most exceptional orbits of the system, and 
not necessarily the most probable behavior of the Solar 
System over its age. For the understanding of the general 
behavior of the Solar System, the density function for the 
eccentricity and inclination of the planets is a comple- 
mentary information that can be very valuable for the 
analysis of many physical parameters during the evolu- 
tion of the Solar System. It was for example useful for 
the understanding of the capture of Mercury into the 3/2 
spin orbit resonance (Correia and Laskar, 2004), or for the 
analysis of the past climate evolution of Mars (Laskar et 
al, 2004b). Here we have systemized the approach elab- 
orated in PI for all the planets of the Solar System, with 
statistics over 1001 different orbits in order to obtain a 
complete statistical view of the variations of the orbital 
parameters of the Solar System. 

Because of the chaotic behavior of the system, we know 
that it will never be possible to retrieve precisely the past 
(or future) orbital evolution of the Solar System over more 
than a few tens of millions of years (Laskar, 1989, 1990). 
On the other hand, beyond about 250 Myr one can obtain 
very smooth density functions for the eccentricity of the 
planets that will tell us the main general behavior of the 
orbits. Moreover, these density functions are some of the 
only accurate informations one can obtain beyond a few 
100 Myr. 

As in PI, we have divided here the time interval in 250 
Myr intervals, and statistics are done over each 250 Myr 
interval, using the set of 1001 orbital solutions in nega- 
tive time with the initial conditions of Table [T] for which 
the output has been recorded with a 1000 yr step size. 
The first 250 Myr interval is discarded, as the random- 
ization due to the chaotic evolution of the system has not 
yet taken place (see PI for a more complete analysis of 
Mars eccentricity solution), and the normalized probabil- 
ity distributions functions (PDF), for the eccentricities 
and inclinations of all the planets are displayed in figures 
0and Hi 

4.1 Discussion 

The difference in the density functions for the eccentric- 
ity (or inclination ) of the inner (Mercury, Venus, Earth, 
Mars) or outer planets (Jupiter, Saturn, Uranus, Nep- 
tune) is striking on these plots. First of all, the shape of 
the PDF is different. For the inner planets, the PDF is 
similar to a Gaussian curve, although slightly different. 
In particular, all the curves have zero values for e = 0, 
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Figure 9: Normalized density functions for the eccentricities of the planets. The statistic is established with 1001 solutions with very close 
initial conditions. The evolution is followed over 19 intervals of 250 Myr, represented by the different curves. The variation of these curves 
thus reflects the chaotic diffusion of the solutions. For the outer planets, all 19 curves are practically identical. 



with a positive, linear slope, while they have some long 
tail for large values of the eccentricity. 

On the opposite, for the outer planets the density curves 
are strictly confined between two non zero values, and 
the curves have two peaks close two the minimum and 
maximum value of the eccentricity (resp. inclination). 

Another striking feature is the fact that for the inner 
planets, we see clearly that numerous curves are displayed 
in each plot, as the density curve is different for each time 
interval of 250 Myr. This is the result of a significant 
diffusion that occurs in the eccentricity and inclination 
(Laskar, 1994), while for the outer, as the diffusion is 
practically non existent, all the 19 density curves that are 
actually plotted in figures [9] and [10] are virtually identical 
and superpose nearly exactly. 

The difference between the two kinds of PDF will even 
be more striking in the next section, as we will attempt to 
associate to each PDF the density function of a simple dy- 
namical system. For the inner planets, we found that the 
densities were very well modclized by a Rice distribution, 
(Rice, 1945) that is characteristic of either the modulus 
of a random walk in two dimensions, or alternatively to 
a quasiperiodic signal with noise (Rice 1945). We will 
denotes these density functions Rice densities. 



On the opposite, the outer planets density function is in 
fact representative of a density function of a quasiperiodic 
signal with a moderate number of periodic terms. We will 
call such a density function a quasiperiodic density. 

In the next section, we will examine more closely these 
two different densities. 

5 Probability density functions 
5.1 Rice densities 

The Rice distribution is a continuous probability distri- 
bution with density function 

X I X 2 + m 2 \ I XTTl \ 

U,m(x) = -j exp I —5— I J I —J- I (2) 

where Iq(z) is the modified Bessel function of the first 
kind obtained with its series expansion 

00 / \ 2 " 1 
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Figure 10: Normalized density functions for the inclinations of the planets. The statistic is established with 1001 solutions with very close 
initial conditions. The evolution is followed over 19 intervals of 250 Myr, represented by the different curves. The variation of these curves 
thus reflects the chaotic diffusion of the solutions. For the outer planets, all 19 curves are practically identical. 



This distribution is obtained for example for the mod- 
ulus of z — x + iy where x,y are two Gaussian indepen- 
dent variables with variance a 2 and mean a. r ,a„, with 



5.2 Diffusion over 5 Gyr 



m - 

These PDF could thus be well suited for the eccen- 
tricity, if the rectangular variables h — esin(tu) and 
k = ecos(zu) become practically Gaussian random vari- 
ables because of the chaotic diffusion. 

A Rice PDF has been successfully fitted to the eccen- 
tricity and inclination of the inner planets for each of 
the 250 Myr time intervals. It is impossible to display 
here graphics demonstrating the quality of the adjustment 
for all of these intervals, so we will reproduce here only 
the most representative cases, for the eccentricity of Mer- 



cury (Figjllj, Venus (Fig|12|, Earth (Fig|13|, and Mars 
(Fig 14 1, at 3 different epochs, while all cases can be sum- 
marized in Table [6] 

It is indeed remarkable that the PDF of the eccentricity 
and inclination of the planets are such smooth functions. 
It is as well remarkable that they are so well approximated 
by simple 2 parameters PDF (f a ,m) of a single curve fam- 
ily with relatively simple expression ([2]). 
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m 


bo 


b x 
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0.1875 


2.070e-03 


1.043e-03 
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0.02235 + 0.00014 T 


4.197e-04 


5.110e-06 


3 


0.01951 + 0.00013T 


3.181e-04 


3.323e-06 


4 


0.06437- 0.00188 T 


1.002e-03 


9.127e-05 


1 


4.9896 


9.040e+00 


1.272e+00 


2 


1.5864 


1.492c+00 


2.377e-02 


3 


1.5803 


1.063e+00 


2.308e-02 


4 


4.8289 - 0. 1703 T 


3.623e+00 


2.398e-01 



Table 6: Evolution of the parameters for the Rice PDF of the inner 
planets Mercury, Venus, the Earth, and Mars over 5 Gyr. N is the 
index of the planet. The parameters m of the Rice PDF (Eqs. [2| is 
given in column 2, while the coefficients t>0i b\ allow one to compute 
the a parameter of (TJjl as a 2 = bo + b\T, where T is in Gyrs. 

As the PDF of eccentricity and inclination are well ap- 
proximated for various epochs by Rice PDF, we have de- 
rived by least square fit, the values of the parameters m, tr, 
of these PDF for all variables and all time intervals of 250 
Myr. Except for Mars, the m value (related to the mean) 
do not present large variations over time (Figs. 15 16 17 
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Figure 11: PDF of the eccentricity of Mercury (full line) for three 
different dates from J2000 (500 Myr, 2500 Myr, 5000 Myr). For 
each case, a two parameters Rice PDF (dashed line) has been fitted 
to the eccentricity PDF. 



Figure 12: PDF of the eccentricity of Venus (full line) for three 
different dates from J2000 (500 Myr, 2500 Myr, 5000 Myr). For 
each case, a two parameters Rice PDF (dashed line) has been fitted 
to the eccentricity PDF. 



18 1. On the opposite, the standard deviation parameter a 
increases with time, and we have a nearly linear relation 



cr 2 = b + bxT 



(4) 



similar to a diffusion process. In Table [6] are gathered the 
fitted values of m, a 1 over 5 Gyr. 

We have thus in equation [2] and Table [6] some very sim- 
ple formulas that will allow one to represent the PDF of 
the eccentricity and inclination of the inner planets over 
very long times, of several Gyr. Quite remarkably, al- 
though it is impossible to predict the precise evolution of 
the individual trajectories, we are able to give very sim- 
ple expressions that fit well with the observed eccentricity 
and inclination PDF of the inner planets. 

These formulas can then be used for the analysis of 
the past or future behavior of the Solar System. They 
could be used for example to compute an estimate of the 
capture probability of Mercury in the 3/2 resonance with- 



out requiring heavy numerical computations (Correia and 
Laskar, 2004). 

We have used here a different PDF than in (Laskar tt 
al. , 2004b) , where a model of a random walk with an ab- 
sorbing edge at zero was used for Mars eccentricity. In 
fact, the results obtained with the two PDF are very sim- 
ilar in this case, but we prefer the Rice formulation that 
can be more easily interpreted. Moreover, we see here that 
the evolution of the system is well describe by a diffusive 
process for the rectangular coordinates k, h that behave 
like random Gaussian variables on long time scale. Actu- 
ally, the variance of the eccentricity evolves linearly with 
time. This is somewhat different from the observation 
made in (Laskar, 2004) for Mars eccentricity, but this is 
because here we are searching for models that fit more 
generally with the behavior of all the inner planets, for 
both eccentricity and inclination. 
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Figure 13: PDF of the eccentricity of the Earth (full line) for three 
different dates from J2000 (500 Myr, 2500 Myr, 5000 Myr). For 
each case, a two parameters Rice PDF (dashed line) has been fitted 
to the eccentricity PDF. 




Figure 14: PDF of the eccentricity of Mars (full line) for three 
different dates from J2000 (500 Myr, 2500 Myr, 5000 Myr). For 
each case, a two parameters Rice PDF (dashed line) has been fitted 
to the eccentricity PDF. 



6 Outer planets 

The PDF of the outer planets (Figs. |9j [lOj are very dif- 
ferent from the PDF of the inner planets. Moreover, as 
it was already stated, the diffusion is practically incxis- 
tent, and all PDF curves over the different time intervals 
are virtually identical. The shape of these outer planets 
PDF can be understood as the PDF of some quasiperi- 
odic functions of time with a small number of harmonics. 
Indeed, for a simple periodic function, g(t) = sin(i), the 
PDF is 



(5) 



where l]_i ^(x) is the characteristic functiorj^] of the in- 
terval ] -1, 1[ (FigMl. 



2 The characteristic function l] a u(t) of the interval ]a, b[ is de- 
fined as 1] a. ,b[ (*) = 1 for t £]a, 6[, and l] a ^(t) = if t £]a, b[. 



One can now understand that when for g(t), instead 
of a single sine term, we have several periodic terms, the 
PDF of g(t) will become slightly distorted with respect 
to a pure sine function, and we will recover the specific 
form of the PDF of the outer planets eccentricities and 
inclinations (Figs. [9 10). In order to illustrate this, we 
will approximate the eccentricity of Jupiter with a few 
periodic terms, obtained through frequency map analysis 
(Laskar, 1990, 2005). 

With a frequency analysis of the eccentricity of Jupiter 
over 50 Myr, we obtain a quasiperiodic approximation of 
the eccentricity with 5 periodic terms on the form 



e = a + di cos(vit + fa) 



(G) 



are given in table [7] The ec- 
centricity of Jupiter over 50 Myr is plotted in figure [20] as 



where the values of a* , Vi , 
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Figure 15: Evolution with time of the parameters m (top) and b 
(bottom) of the PDF of the eccentricity of Mercury, m is approx- 
immated by a constant value mo (dotted line), while b is fitted with 
a linear slope b = bo + biT, where the time T is in Gyr. The fitted 
values mo , bo , bi are given in Table [6] 
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Figure 16: Evolution with time of the parameters m (top) and b 
(bottom) of the PDF of the eccentricity of Venus, m and b are fitted 
with a linear slope m = mo + m\T and b = bo + feiT respectively, 
where the time T is in Gyr. The fitted values mo, mi, bo, fei are 
given m Table [6] 
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Table 7: Quasiperiodic decomposition of the eccentricity of Jupiter 
over 50 Myr. we have e = ao + ^2 a; cos^t + (f>i). The decompo- 
sition is obtained through frequency analysis (Laskar, 1988, 2005). 
The residuals are given in figure 1201 



well as its difference with this quasiperiodic approxima- 
tion. The density function obtained with this quasiperi- 



odic approximation is plotted in figure 21 in dotted line 
together with the density function of the full numerical 
integration (in full line). The two PDF are very close. It 
will be the same for the other outer planets. For these 
planets with a motion that is very close to quasiperiodic, 
the quasiperiodic decomposition obtained over 50 Myr as 
in (Laskar, 1990) provides in fact a very good representa- 
tion of the orbit over several Gyr. It is thus not necessary 
to reproduce these expressions here, and we would rather 
provide in a forthcoming publication some full and accu- 
rate quasiperiodic representations for the motion of the 
outer planets, in agreement with our latest adjustment 
of the planetary ephemeris to observations (Fienga et at, 
2007). 



7 Conclusions 

In (Laskar 1994), I showed the possibility of a very large 
increase of the eccentricity of Mercury, allowing for a close 
encounter or a collision with Venus. But in this work, 
there was no estimate of the probability for such an event 
to take place within a few Gyr. Here, such an estimate 
is given, by the extensive study of more than 1001 orbits. 
Indeed, it is found that the probability for Mercury's ec- 
centricity to exceed 0.6 within 5 Gyr is about 1 to 2 % 
which can be considered as a large value for such an im- 
portant event. 

When the contributions of general relativity and the 
Moon are not taken into account, this probability in- 
creases in a large amount and in the present numerical 
simulations, nearly half of the orbits went beyond 0.6 in 
less than 4 Gyr. As this appear to be in contradiction 
with the numerical results of (Ito and Tanikawa 2002), I 
have also performed a direct numerical simulation of the 
Newtonian equations, with 10 nearby initial conditions, 
without the Moon or general relativity over 3 Gyr, and 
found results that are in good agreement with the results 
of the integration of the secular system. In particular, 
two orbit were found for which the eccentricity of Mer- 
cury rises to very large values, beyond 0.8, thus allowing 
for a close encounter or a collision with Venus. 

This direct numerical simulation is performed with 
slightly better accuracy than the numerical integration of 
(Ito and Tanikawa 2002). One should wonder about the 
reason of such a different behavior in the present compu- 
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Figure 17: Evolution with time of the parameters m (top) and 
6 (bottom) of the PDF of the eccentricity of the Earth, m and 
b are fitted with a linear slope m = mo + m\T and b = bo + 
b\T respectively, where the time T is in Gyr. The fitted values 
mo, mi, bo , fei are given in Table |6] 
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Figure 18: Evolution with time of the parameters m (top) and b 
(bottom) of the PDF of the eccentricity of Mars, m and b are fitted 
with a linear slope m = mo + m\T and b = bo + 6iT respectively, 
where the time T is in Gyr. The fitted values mo,mi,6oi^i are 
given in Table |6] 



tations and in the simulation of (Ito and Tanikawa 2002) . 
The most probable reason for this is that the two inte- 
grations do not have the same initial conditions or model. 
The solutions of (Ito and Tanikawa 2002) may evolve in 
a slightly more regular region of the phase space. On the 
other hand, the secular equations are in very good agree- 
ment with the present direct integrations. This comfort 
their reliability and usefulness as for the secular equations, 
it was possible to perform an extended statistical study 
over 1001 solutions. One could perform additional direct 
numerical integrations in order to verify the probability 
law for the excursion of Mercury's eccentricity, but this is 
of limited interest if it concerns the model of pure Newto- 
nian equations, and not the full model including general 
relativity. 

Indeed, as we have demonstrated here in section |3.4| 
the contribution of general relativity changes in a consid- 
erable manner the behavior of the Solar System dynam- 
ics. Indeed, in absence of GR, the secular frequency of 
the perihelion of Mercury gi decrease by 0.43"/yr, and 
becomes closer to the secular frequency of the perihelion 
of Jupiter (75. As the eccentricity varies under the secu- 
lar planetary perturbations, the system can enter into the 
<?i — .95 secular resonance that can drive Mercury's eccen- 
tricity to very high values, beyond 0.8, where additional 
short period chaotic behavior occurs, inducing changes of 
semi-major axis of the planet (Fig. |8|. This is indeed the 
mechanism that is present in our two solutions of the New- 
tonian equations 52,5s for which Mercury's eccentricity 
increased beyond 0.8 (Fig. [7]). 



An important result of the present paper is to show this 
possibility of very large increase of Mercury's eccentricity, 
beyond 0.8, allowing for a possible collision with Venus. It 
was actually possible to see that with such a large value 
of the eccentricity, the planet entered a region of mean 
motion resonances and chaos inducing changes in its semi- 
major axis. 

It is clearly desirable to conduct now a full scale nu- 
merical experiment with the complete equations includ- 
ing general relativity and the contribution of the Moon, 
in order to provide some precise results on the chaotic 
behavior and possible evolution of our Solar System. We 
are indeed planning such a numerical study for the near 
future. 

Apart from demonstrating the possibility of very large 
values for the eccentricity of Mercury, as was forecasted 
in (Laskar, 1994), I provide here probability density func- 
tions (PDF) for the eccentricities and inclinations of the 
terrestrial planets. It is in fact quite remarkable that the 
PDF of all terrestrial planets, when computed over long 
time, beyond about 500 Myr, become very smooth func- 
tions that are well approximated by a very simple two 
parameters function, namely the Rice distribution. This 
provides very compact formulas for the evolution of these 
PDF over time for all the terrestrial planets. The evo- 
lution of the PDF with time is similar to the one of a 
diffusion process with a linear increase of the variance 
variable a 2 with time. 

For the outer planets, we have seen here that the chaotic 
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Figure 19: Plot of the function f(x) = l[_i,i] (x)/(n\/l — x 2 ) that 
is the PDF of g(t) = sin(t) (Eqj5|. 




Figure 20: Solution for the eccentricity of Jupiter (top), and resid- 
uals after substraction of a quasiperiodic approximation with 5 pe- 
riodic terms (Eq. |Bj. 



diffusion is extremely small, without practical change of 
the PDF over time. The PDF being well approximated 
by the PDF of a quasiperiodic approximation of the so- 
lution. In the full equations of motion of the outer So- 
lar System, there is some chaotic behavior that has been 
described (Sussman and Wisdom, 1992, Murray and Hol- 
man, 1999, Guzzo, 2005), but although a more detailed 
analysis should be made, It seems that the diffusion in- 
duced by this intrinsic chaotic behavior resulting from 
mean motion interactions should be smaller than the dif- 
fusion induced by the chaotic behavior of the inner Solar 
System that is already included in the present work. It 
is thus doubtful that the PDF for the outer planets that 
have been computed here will be much changed by the 
consideration of the full dynamics of the outer Solar Sys- 
tem. 

More generally, the PDF that are given here are the 
PDF obtained with the secular equations. Nevertheless, 
they should be very stable with respect to small changes 
in the model or in the involved parameters, and we can 
conjecture as well that they are very close to the PDF of 
the full system. 




0.04 0.06 
eccentricity 



Figure 21: Density function for the eccentricity of Jupiter (solid 
line) and in dotted line, the density function of a quasiperiodic ap- 
proximation with 5 periodic terms (Eq. pi. 
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